library(hydroGOF)
library(clusterSim)
file = "myd"
data=read.csv(paste("C:/in/",file,"_merge.csv",sep=""))
data=subset(data,data$station!=3)

data$aod=data.Normalization(data$aod,type="n5",normalization="column")
data$temp=data.Normalization(data$temp,type="n5",normalization="column")
data$avg_temp=data.Normalization(data$avg_temp,type="n5",normalization="column")
data$avg_rh=data.Normalization(data$avg_rh,type="n5",normalization="column")
data$avg_preci_24=data.Normalization(data$avg_preci_24,type="n5",normalization="column")


model=lm(pm25~aod+temp+avg_temp+avg_rh+avg_preci_24,data)
nSample=nrow(data)
nPredict= length(model$coefficients)-1

cutoff=4/(nSample-nPredict-1)
cook_thresold=qf(0.2,nPredict+1,nSample-nPredict-1)

cookValue=cooks.distance(model)
data["cookValue"] = cookValue


data1=data
data2=subset(data,data$cookValue<=cutoff)

model22=lm(pm25~aod+temp+avg_temp+avg_rh+avg_preci_24,data2)
nSample22=nrow(data2)
nPredict22= length(model22$coefficients)-1

cutoff2=4/(nSample22-nPredict22-1)
cookValue22=cooks.distance(model22)
data2["cookValue"] = cookValue22

data3=subset(data,data$cookValue<=cook_thresold)

result=data.frame()

for (i in 1:5) {
	test_station = i
	trainData1=subset(data1,data1$station!=test_station)
	testData1=subset(data1,data1$station==test_station)
		
	model1=lm(pm25~aod+temp+avg_temp+avg_rh+avg_preci_24,trainData1)
	pm25model=predict(model1,testData1)
	
	model_re1=sum(abs(testData1$pm25-pm25model)/testData1$pm25)/nrow(testData1)*100
	model_rmse1=rmse(testData1$pm25,pm25model)
	model_r1=cor(testData1$pm25,pm25model)
	
	total_train1=nrow(trainData1)
	total_test1=nrow(testData1)
	r1=round(model_r1*model_r1,3)
	rmse1=round(model_rmse1,3)
	re1=round(model_re1,3)
	
	print(paste("Year model 1: ",test_station))
	print(paste("Number of training: ",total_train1))
	print(paste("Number of testing: ",total_test1))
	print(paste("R2:",r1))
	print(paste("RMSE:",rmse1))
	print(paste("RE:",re1))
	
		
	print("------------------")
	
	
	trainData2=subset(data2,data2$station!=test_station)
	testData2=subset(data2,data2$station==test_station)
	
	model2=lm(pm25~aod+temp+avg_temp+avg_rh+avg_preci_24,trainData2)
	pm25model=predict(model2,testData2)
	
	
	model_re2=sum(abs(testData2$pm25-pm25model)/testData2$pm25)/nrow(testData2)*100
	model_rmse2=rmse(testData2$pm25,pm25model)
	model_r2=cor(testData2$pm25,pm25model)
	
	total_train2=nrow(trainData2)
	total_test2=nrow(testData2)
	r2=round(model_r2*model_r2,3)
	rmse2=round(model_rmse2,3)
	re2=round(model_re2,3)
	
	print(paste("Year model 2: ",test_station))
	print(paste("Number of training: ",total_train2))
	print(paste("Number of testing: ",total_test2))
	print(paste("R2:",r2))
	print(paste("RMSE:",rmse2))
	print(paste("RE:",re2))
	
	print("------------------")
	

		
	trainData3=subset(data3,data3$station!=test_station)
	testData3=subset(data3,data3$station==test_station)
	
	model3=lm(pm25~aod+temp+avg_temp+avg_rh+avg_preci_24,trainData3)
	pm25model=predict(model3,testData3)
	
	
	model_re3=sum(abs(testData3$pm25-pm25model)/testData3$pm25)/nrow(testData3)*100
	model_rmse3=rmse(testData3$pm25,pm25model)
	model_r3=cor(testData3$pm25,pm25model)
	
	total_train3=nrow(trainData3)
	total_test3=nrow(testData3)
	r3=round(model_r3*model_r3,3)
	rmse3=round(model_rmse3,3)
	re3=round(model_re3,3)
	
	print(paste("Year model 2: ",test_station))
	print(paste("Number of training: ",total_train3))
	print(paste("Number of testing: ",total_test3))
	print(paste("R2:",r3))
	print(paste("RMSE:",rmse3))
	print(paste("RE:",re3))
	
	print("------------------")
	result=rbind(result,data.frame(test_station,total_train1,total_test1,r1,rmse1,re1,total_train2,total_test2,r2,rmse2,re2,total_train3,total_test3,r3,rmse3,re3))

}
write.csv(result,paste("C:/a/",file,"_station.csv",sep=""))

